Recap

  • What is cluster analysis?
  • Distance Measure
  • K-means algorithm
  • Hierarchical algorithms
  • Dendrograms

Outline

  1. What is Modeling?
  2. Correlation
  3. Simple Linear Regression
  4. Classical Normal Linear Regression Model
  5. Regression Model Diagnostics
  6. Model Selection

Modeling

Working with A Single Variable

When analyzing a single numerical variable, we can use numerical statistics to gain insights into key attributes such as

  • central tendency, e.g. mean, median, and mode
  • dispersion, e.g. variance, standard deviation and IQR
  • shape, e.g. skewness and kurtosis

Alternatively, we can use graphical tools like box plots and density plots to represent these characteristics.

maxt <- read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%
  pull(maxt)

p1 <- ggplot() +
  geom_boxplot(aes(maxt)) +
  xlab(expression(degree*C)) +
  theme_light()

p2 <- ggplot() +
  geom_density(aes(maxt)) +
  xlab(expression(degree*C)) +
  theme_light()

patchwork::wrap_plots(p1, p2) +
  patchwork::plot_annotation(title = "Maximum daily temperature for bushfires in Victoria between 2015 and 2018")

Bivariate Relationships

When working with two numerical variables, a scatterplots is a natural choice to explore association between them.

It helps identify whether the variables are:

  • Positively associated: as one variable increases, the other also increases
  • Negatively associated: as one variable increases, the other decreases
  • No association: the data points are scattered randomly, showing no clear relationship

read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%
  ggplot() +
  geom_point(aes(mint, maxt), alpha = 0.6) +
  xlab(expression(Minimum~temperature~(degree*C))) +
  ylab(expression(Maximum~temperature~(degree*C))) +
  theme_light() +
  ggtitle("Daily maximum and minimum temperature for bushfires \nin Victoria between 2015 and 2018")

Correlation

Correlation is the linear association between two variables, ranging from \(-1\) to \(1\).

For two variables, \(X\) and \(Y\), the correlation coefficient \(r\) is calculated as follows:

\[r_{xy} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n}(y_i - \bar{y})^2}} = \frac{S_{xy}}{S_x S_y},\]

where \(S_{xy}\) represents the biased sample covariance between \(X\) and \(Y\), and \(S_x\) and \(S_y\) denote the biased standard deviations of \(X\) and \(Y\), respectively.

Tip

  • In correlation, \(x\) and \(y\) are interchangeable, i.e. \(r_{xy} = r_{yx}\).

Correlation

Interpretation:

  • \(r = 1\): strong positive correlation
  • \(r = -1\): strong negative correlation
  • \(r = 0\): variables are not linearly associated, though this doesn’t imply they have no association!

map_df(seq(1, -1, -0.2), function(r) {
  mvtnorm::rmvnorm(500, sigma = matrix(c(1, r, r, 1), ncol = 2)) %>%
    as.data.frame() %>%
    mutate(r = r)
}) %>%
  ggplot() +
  geom_point(aes(V1, V2)) +
  facet_wrap(~r) +
  theme_light() +
  xlab("X") +
  ylab("Y")

Correlation vs Causation

Causation means that one variable directly influences another.

Here are some examples:

  • Smoking increases the risk of lung cancer.
  • Drinking alcohol impairs motor skills and cognitive function.

Remember: Correlation does not imply causation!

  • Both ice cream sales and drowning incidents rise in hot weather, but ice cream sales do not cause drownings. The real factor is the temperature.
  • As shoe size increases, so does a child’s reading ability. However, this relationship is actually due to age.
  • The number of pirates and global temperatures show a correlation, but this is a spurious relationship — two unrelated variables that happen to coincide!

What is Modeling?

Modeling is a method for understanding the relationship between different variables.

We typically represent the relationship of one response variable, \(Y\), in relation to other explanatory variables, \(X_1, X_2, \ldots, X_p\), using a mathematical function:


\[Y = f(X_1, X_2, ..., X_p).\]


Why we need modeling?

  • While we can observe the data, we may not have insights into the underlying processes, i.e. the joint distributions of the variables.

  • Therefore, we rely on data to make inferences about these processes and build models to explain them.

Different Types of Models

There are lots of different models.

  • Linear model
  • Non-linear model
  • Decision tree
  • Neural network
  • Non-parametric model

For now we focus on linear models.

read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%
  ggplot() +
  geom_point(aes(mint, maxt)) +
  geom_smooth(aes(mint, maxt), method = "lm", se = FALSE) +
  xlab(expression(Minimum~temperature~(degree*C))) +
  ylab(expression(Maximum~temperature~(degree*C))) +
  theme_light() +
  ggtitle("Daily maximum and minimum temperature for bushfires \nin Victoria between 2015 and 2018")

Simple Linear Regression

A simple linear regression (SLR) is a regression model that involves a single explanatory variable. The formula is expressed as follows:

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad i = 1, \ldots, n.\]

  • \(x_i\) and \(y_i\): the \(i\)-th observation’s \(x\) and \(y\) value
  • \(\beta_0\) and \(\beta_1\): unknown constants that represent the intercept and slope, also known as coefficients or parameters.
  • \(\varepsilon\): the error term, which has specific properties and this is the stochastic component of the model, enabling us to make statistical inferences about the model and its parameters!

Visulize Simple Linear Regression with ggplot

What are the intercept and slope?

Visulize Simple Linear Regression with ggplot

What are the intercept and slope?

Fitting Simple Linear Regression with lm

\[\text{maxt}_i = \beta_0 + \beta_1\text{mint}_i + \varepsilon_i, \quad i = 1, \ldots, n.\]

  • \(\hat{\beta}_0 = 12.02\) and \(\hat{\beta}_1 = 1.04\) are the estimated values of the coefficients.

Fitting Simple Linear Regression with lm

\[\text{speed}_i = \beta_0 + \beta_1\text{dist}_i + \varepsilon_i, \quad i = 1, \ldots, n.\]

  • \(\hat{\beta}_0 = -17.58\) and \(\hat{\beta}_1 = 3.93\) are the estimated values of the coefficients.

Estimation

But how are these coefficients actually estimated?

First, we need to introduce the concepts of fitted values and residuals.

\[\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i.\]

  • The hat symbol ( \(\hat{}\) ) indicates that these values are estimated or fitted.
  • \(\hat{\beta}_0\) and \(\hat{\beta}_1\) are the estimated coefficients.
  • \(\hat{y}_i\) is the \(i\)-th fitted value.

The \(i\)-th residual \(e_i\) is the difference between the fitted value and the actual value of the response variable:

\[ e_i = y_i - \hat{y}_i = y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i). \]

Estimation

Estimation

We aim to minimize the residuals, meaning we want the regression line to be as close to the data points as possible.

But how do we determine which solution is better? How do we evaluate it?

Estimation

We can minimize the mean squared error (MSE) by solving for the coefficients:

\[ \underset{\beta_0, \beta_1}{\text{argmin}} \, \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \]

This ensures that the regression line would be the line closest to the cloud of points!

Estimation

The minimization process described before leads to the following estimates:

\[\hat{\beta}_1 = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i -\bar{x})^2} \Rightarrow \text{slope}, \text{ and}\]

\[\hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x} \Rightarrow \text{intercept},\] where \(\bar{y} = \frac{1}{n}\sum_{i=1}^ny_i\) and \(\bar{x} = \frac{1}{n}\sum_{i=1}^nx_i\) are the sample mean.

Note

Recall that \[r_{xy} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n}(y_i - \bar{y})^2}}.\]

So, for simple linear regression, \(\hat{\beta}_1 = r_{xy}\frac{S_y}{S_x}\).

Prediction

Once the model is fitted and all parameters are estimated, we can make predictions by plugging in \(x\) to obtain the predicted value \(\hat{y}\).

This predicted value reflects the average \(y\) for each corresponding value of \(x\), i.e. \(\hat{E}[y|x]\).

\[\widehat{\text{maxt}}_i = 12.02 + 1.04 \times \text{mint}_i.\]

Examples:

  • For \(\text{mint} = 10\), \(\widehat{\text{maxt}} = 10 \times 1.04 + 12.02 = 22.42\). So, on average, we expect a day with a minimum temperature of 10 to have a maximum temperature of 22.42.

Warning: While we “expect” this to happen, there will be some variability, which is captured by the error term.

  • Predictions are generally valid for new values within the range of the original data used to fit the model. Extrapolation beyond this range is less reliable and carries greater risk.

Coefficient Interpretation

Note that, when \(\text{mint} = 0\), the predicted \(\widehat{\text{maxt}} = 0 \times 1.04 + 12.02 = 12.02\).

  • This value is equal to the intercept, indicating that the intercept represents the predicted average maximum temperature when the minimum temperature is zero.

Similarly, the slope \(\hat{\beta}_1 = 1.04\) can be interpreted as follows:

  • For each additional 1 degree increase in minimum temperature, the maximum temperature is expected to increase, on average, by 1.04 degrees.

Classical Normal Linear Regression Model

We have discussed the most basic regression model in statistics. Now, let’s consider its extension to include multiple explanatory variables:

\[y_i = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... +\beta_px_p + \varepsilon_i, \quad i=1,...n,\]

where \(\varepsilon_i \sim N(0, \sigma^2)\) for \(i = 1,...n\).

Model assumptions:

  1. Linearity: The relationship between explanatory variables and response variable is linear.
  2. Independence: The errors are assumed to be independent of each other.
  3. Homoscedasticity: The errors have constant variance across all levels of the explanatory variables.
  4. Normality: The errors are assumed to be normally distributed.

Classical Normal Linear Regression Model

Coefficient Interpretation

In a CNLRM, the interpretation of the intercept term \(\hat{\beta}_0\) remains the same.

  • It represents the predicted average value of \(y\) when all \(x\) values are equal to zero.

For the other parameters \(\hat{\beta}_j\) where \(j = 1, \ldots, p\), the interpretation is as follows:

  • While holding all other variables constant, for each additional unit increase in \(x_j\), the variable \(y\) is expected to increase, on average, by \(\hat{\beta}_j\) units.

Let’s have a break!

Model Diagnostics

Model Diagnostics

All models are wrong; some are useful - George Box.

After fitting a linear model, it is crucial to assess whether the model assumptions are satisfied.

There is a real risk, that a model is imposing structure that is not really there.

It can result in inaccurate estimates of the coefficients and misleading interpretations of the relationship between \(x\) and \(y\).

Residual plots

The first thing to check is the residuals versus fitted values plot.

What constitutes a good residual plot?

  • Residuals should be randomly scattered around the horizontal line at zero.
  • There should be no discernible patterns in the residuals, such as curves or trends.
  • The spread of the residuals should remain consistent across all levels of the fitted values.

Residual versus Explanatory Variable

The second thing to check is the residuals versus each individual explanatory variable plot.

The principles for a good residual plot remain the same.

This type of plot offers clearer insights into the relationship between the residuals and each individual explanatory variable, which can be invaluable for refining and improving the model.

Q-Q plot

The third thing to check is the Q-Q plot.

To check a Q-Q plot, compare the distribution of your residuals to a theoretical normal distribution.

In a good Q-Q plot, the points should closely follow the 45-degree reference line.

Significant deviations from this line, such as curves or outliers, indicate that the residuals may not be normally distributed, which could suggest violations of model assumptions.

Lineup of Residual Plots

We can display the actual residual plot alongside a matrix of residual plots simulated under the assumption that the model is correct.

If the true residual plot is distinct and identifiable within the matrix, it provides evidence suggesting a violation of the model assumptions.

This method is known as the lineup protocol.

  • Learn more about it from here.

# remotes::install_github("autovi")
mod <- lm(maxt ~ mint + se, data = bushfire)
autovi::auto_vi(mod)$plot_lineup(alpha = 0.1)

Model Selection

Beyond checking model assumptions, we often need to select between different models. For instance, choosing between maxt ~ rf + mint or maxt ~ mint.

There are several statistics to assess model fit:

  • \(R^2\)
  • adjusted \(R^2\)
  • AIC
  • BIC
  • Deviance

\(R^2\)

\[R^2 = 1-\frac{\sum_{i=1}^{n}(y_i-\hat{y})^2}{\sum_{i=1}^{n}(y_i-\bar{y})^2} = \frac{\text{Residual sum of square}}{\text{Total sum of square}} \in [0,1].\]

  • The coefficient of determination, \(R^2\), ranges from 0 to 1, with 1 indicating a perfect fit.
  • Adding more variables increases \(R^2\) but also adds model complexity.
  • Adjusted \(R^2\) penalizes for excessive variables to account for model complexity.
  • \(R^2\) is a common measure of linear model fit, indicating the percentage of variability in \(y\) explained by the model, with the remainder attributed to factors not included.

\(R^2\)

Which model is better?

AIC, BIC, and Deviance

AIC, BIC, and Deviance are goodness of fit measure to compare models.

  • AIC = Akaike Information Criterion (can be used to compare models. The smaller the value the better the model.)

  • Similarly BIC = Bayes Information Criterion indicates how well the model fits, best used to compare two models. Lower is better.

  • Deviance is the residual variation, how much variation in response that IS NOT explained by the model. The close to 0 the better, but it is not on a standard scale. In comparing two models if one has substantially lower deviance, then it is a better model.

Going beyond A Single Model

Gapminder

  • Hans Rosling was a Swedish doctor, academic and statistician, Professor of International Health at Karolinska Institute. Sadly he passed away in 2017.

  • He developed a keen interest in health and wealth across the globe, and the relationship with other factors like agriculture, education, energy.

  • You can play with the gapminder data using animations at https://www.gapminder.org/tools/.

Watch the Video

R package: gapminder

Contains subset of the data on five year intervals from 1952 to 2007.

Change in Life Expectancy in Countries over Time

Change in Life Expectancy in Countries over Time

  • There generally appears to be an increase in life expectancy
  • A number of countries have big dips from the 70s through 90s
  • A cluster of countries starts off with low life expectancy but ends up close to the highest by the end of the period.

Fit Linear Models for Multiple Countries with nest

Differences in Life Expectancy by Continent

Fit An Overall Model

  • Categorical explanatory variables are typically encoded to dummy variables, one for each of the levels.
  • Each coefficient describes the expected difference compared to the baseline level.